function [E,EE,legend,condn] = collect_data(output_file,indexes,varargin)
%[E,EE,legend,condn] = collect_data(output_file,indexes,ndigits,expnot)
%
% Collect data from a series of files.
%  output_file: name of the output files using '#' for the numbering
%  indexes: array with the values to be used in output_file as '#'
%  ndigits (optional): number of digits in the file names (default 3)
%  expnot (optional): use 1 to used undexes in exponential notation
%
%  E: struct with the errors
%  EE: same as E but arranged in a single vector; one can use
%        (log(EE(:,1:end-1))-log(EE(:,2:end)))./log(indexes(2:end)./indexes(1:end-1))
%      to compute the convergence rates
%  legend: legend associated to the rows of EE
%  condn: conditioning of the linear system

 if(nargin==2)
   ff = '%.3i'; % default three digits
 elseif(nargin==3)
   ff = ['%.',num2str(varargin{1}),'i'];
 elseif(nargin==4)
   ff = ['%.',num2str(varargin{1}),'e'];
 end

 for i=1:length(indexes)

   istring = num2str(indexes(i),ff);
   pos = strfind(output_file,'#');
   fname = [output_file(1:pos-1),istring,output_file(pos+1:end)];

   % read the file
   vars = load(fname);

   if(isfield(vars,'nnz_mmm11')) E.nnz_mmm11(i)=vars.nnz_mmm11; end
   if(isfield(vars,'e_l2s'))     E.l2e    (i)=vars.e_l2s;     end
   if(isfield(vars,'e_l2k'))     E.l2k    (i)=vars.e_l2k;     end
   if(isfield(vars,'e_pl2s'))    E.pl2e   (i)=vars.e_pl2s;    end
   if(isfield(vars,'e_pl2k'))    E.pl2k   (i)=vars.e_pl2k;    end
   if(isfield(vars,'e_u_l2'))    E.u_l2   (i)=vars.e_u_l2;    end
   if(isfield(vars,'e_u_h1'))    E.u_h1   (i)=vars.e_u_h1;    end
   if(isfield(vars,'e_u_div'))   E.u_div  (i)=vars.e_u_div;   end
   if(isfield(vars,'e_p_l2'))    E.p_l2   (i)=vars.e_p_l2;    end
   if(isfield(vars,'e_p_h1'))    E.p_h1   (i)=vars.e_p_h1;    end
   if(isfield(vars,'e_us_l2'))   E.us_l2  (i)=vars.e_us_l2;   end
   if(isfield(vars,'e_usef_l2')) E.usef_l2(i)=vars.e_usef_l2; end
   if(isfield(vars,'e_q_l2'))    E.q_l2   (i)=vars.e_q_l2;    end
   if(isfield(vars,'e_q_hdiv'))  E.q_hdiv (i)=vars.e_q_hdiv;  end
   if(isfield(vars,'e_qsdg_l2'))  E.qsdg_l2  (i)=vars.e_qsdg_l2;  end
   if(isfield(vars,'e_qsdg_hdiv'))E.qsdg_hdiv(i)=vars.e_qsdg_hdiv;end
   if(isfield(vars,'e_qsrt_l2'))  E.qsrt_l2  (i)=vars.e_qsrt_l2;  end
   if(isfield(vars,'e_qsrt_hdiv'))E.qsrt_hdiv(i)=vars.e_qsrt_hdiv;end
   if(isfield(vars,'e_flux_l2')) E.flux_l2(i)=vars.e_flux_l2;  end

   if(nargout>3) % compute the conditioning
     if(isfield(vars,'mmm'))
       condn(i) = condest( vars.mmm(vars.dofs_nat+1,vars.dofs_nat+1) );
     end
   end

 end

 % for convenience, we also collect the results into one array
 if(exist('E')) % in case we haven't computed the errors
   i = 0;
   for [val, key] = E
     i = i+1;
     legend{i} = key;
     EE(i,:) = val;
   end
 endif

return

